/* Copyright (c) Colorado School of Mines, 2011.*/
/* All rights reserved.                       */

/* SUKDMDCS: $Revision: 1.3 $ ; $Date: 2011/11/16 17:45:18 $	*/

#include "su.h"
#include "segy.h"
#include "header.h"

/*********************** self documentation **********************/
char *sdoc[] = {
" 									",
"  SUKDMDCS - 2.5D datuming of sources for prestack common receiver 	",
" 	     data, using constant-background data-mapping formula.      ",
"            (See selfdoc for specific survey geometry requirements.)   ",
" 									",
"    sukdmdcs  infile=  outfile=  [parameters] 		         	",
"									",
" Required parameters:							",
" infile=stdin		file for input seismic traces			",
" outfile=stdout	file for output  	                        ",
" ttfile		file for input traveltime tables		",
"                                                                       ",
" Required parameters describing the traveltime tables:	         	",
" fzt 			first depth sample in traveltime table		",
" nzt 			number of depth samples in traveltime table	",
" dzt			depth interval in traveltime table		",
" fxt			first lateral sample in traveltime table	",
" nxt			number of lateral samples in traveltime table	",
" dxt			lateral interval in traveltime table		",
" fs 			x-coordinate of first source in table		",
" ns 			number of sources in table			",
" ds 			x-coordinate increment of sources in table	",
"                                                                       ",
" Parameters describing the input data:                                 ",
" nxso                  number of shots                                 ",
" dxso                  shot interval                                   ",
" fxso=0                x-coordinate of first shot                      ",
" nxgo                  number of receiver offsets per shot             ",
" dxgo                  receiver offset interval                        ",
" fxgo=0                first receiver offset                           ",
" dt= or from header (dt)       time sampling interval of input data    ",
" ft= or from header (ft)       first time sample of input data         ",
" dc=0                  flag for previously datumed receivers:          ",
"                          dc=0 receivers on recording surface          ",
"                          dc=1 receivers on datum                      ", 
"                                                                       ",
" Parameters descrbing the domain of the problem:	                ",
" dzo=0.2*dzt		vertical spacing in surface determination	",
" offmax=99999		maximum absolute offset allowed             	",
"                                                                       ",
" Parameters describing the recording and datuming surfaces:            ",
" recsurf=0             recording surface (horizontal=0, topographic=1) ",
" zrec                  defines recording surface when recsurf=0        ",
" recfile=              defines recording surface when recsurf=1        ",
" datsurf=0             datuming surface (horizontal=0, irregular=1)    ",
" zdat                  defines datuming surface when datsurf=0         ",
" datfile=              defines datuming surface when datsurf=1         ",
"                                                                       ",
" Parameters describing the extrapolation:                              ",
" aperx=nxt*dxt/2  	lateral aperture         			",
" v0=1500(m/s)		reference wavespeed             		",
" freq=50               dominant frequency in data, used to determine   ",
"                       the minimum distance below the datum that       ",
"                       the stationary phase calculation is valid.      ",
" scale=1.0             user defined scale factor for output            ",
" jpfile=stderr		job print file name 				",
" mtr=100  		print verbal information at every mtr traces	",
" ntr=100000		maximum number of input traces to be datumed	",
"									",
"                                                                       ",
NULL};
/*
 * Computational Notes:                                                
 *   
 * 1. Input traces must be SU format and organized in common receiver gathers.
 *    
 * 2. Traveltime tables were generated by program rayt2d (or equivalent)     
 *    on any grid, with dimension ns*nxt*nzt. In the extrapolation process,       
 *    traveltimes are interpolated into shot/geophone locations and     
 *    output grids.                                          
 *
 * 3. If the offset value of an input trace is not in the offset array     
 *    of output, the nearest one in the array is chosen.                   
 *
 * 4. Amplitudes are computed using the constant reference wavespeed v0.  
 *                                
 * 5. Input traces must specify source and receiver positions via the header  
 *    fields tr.sx and tr.gx.             
 *
 * 6. Recording and datuming surfaces are defined as follows:  If recording
 *    surface is horizontal, set recsurf=0 (default).  Then, zrec will be a
 *    required parameter, set to the depth of surface.  If the recording  
 *    surface is topographic, then set recsurf=1.  Then, recfile is a required
 *    input file.  The input file recfile should be a single column ascii file
 *    with the depth of the recording surface at every surface location (first 
 *    source to last offset), with spacing equal to dxgo. 
 * 
 *    The same holds for the datuming surface, using datsurf, zdat, and datfile.
 *
 *
 * Assumptions and limitations:
 *
 * 1. This code implements a 2.5D extraplolation operator that allows to
 *    transfer data from one reference surface to another.  The formula used in
 *    this application is an adaptation of Bleistein & Jaramillo's 2.5D data
 *    mapping formula for receiver extrapolation.  This is the result of a
 *    stationary phase analysis of the data mapping equation in the case of
 *    a constant input receiver location (receiver gather). 
 * 
 *
 * Credits:
 * 
 * Authors:  Steven D. Sheaffer (CWP), 11/97 
 *
 *
 * References:  Sheaffer, S., 1999, "2.5D Downward Continuation of the Seismic
 *              Wavefield Using Kirchhoff Data Mapping."  M.Sc. Thesis, 
 *              Dept. of Mathematical & Computer Sciences, 
 *              Colorado School of Mines.
 *
 *
 */
/**************** end self doc ***********************************/
/* Prototype */
  void sum2(int nx,int nz,float a1,float a2,float **t1,float **t2,float **t);
  void dat2d(float *trace,int nt,float ft,float dt,float sx,float gx,
             float **dat,float aperx,int nxso,float fxin,float dxso,float nzo,
             float fzo,float dzo,
             int nxi,float fxi,float dxi,
             float **tsum,int sgn, int nzt,float fzt,float dzt,int nxt,
             float fxt,float dxt,float **szif,float v0,int dc,int ngl,
             float freq);
  void filt(float *trace,int nt,float dt,int sgn);


/* segy trace */
segy tr, tro;

int
main (int argc, char **argv)
{
	int   nt,nzt,nxt,nzo,nxso,nxgo,ns,is,io,ixo,it,
	      ntr,jtr,ktr,mtr,*ng,nxi,sgn,dc,ngl;
	long  nseek;
	float ft,fzt,fxt,fzo,fxso,fxgo,fs,dt,dzt,dxt,dzo,dxso,scale,
	      dxgo,ds,ext,ezt,ezo,exso,exgo,es,s,scal,fxin,as,res;
	float v0,offmax,aperx,sx,gx,
	      dxi,fxi,fgl,freq;
        float ***datg,***ttab,**tsum,**szif;
	
	char  *datain="stdin",*dataout="stdout",*ttfile,*jpfile;
	FILE  *infp,*outfp,*ttfp,*jpfp,*hdrfp;

        int recsurf,datsurf,i;
        float zrec,zdat,r1,r2;

        char  *recfile,*datfile;
        FILE  *recfp=NULL,*datfp=NULL;

    
        /* explicitly define downward continuation */
        sgn=1;


	/* hook up getpar to handle the parameters */
	initargs(argc, argv);
	requestdoc(1);

	/* open input and output files	*/
	if( !getparstring("datain",&datain)) {
		infp = stdin;
	} else  
		if ((infp=fopen(datain,"r"))==NULL)
			err("cannot open datain=%s\n",datain);
	if( !getparstring("dataout",&dataout)) {
		outfp = stdout;
	} else  
		outfp = fopen(dataout,"w");
	fseek(infp,0,1);
	fseek(outfp,0,1);

	if( !getparstring("ttfile",&ttfile))
		err("must specify ttfile!\n");
	if ((ttfp=fopen(ttfile,"r"))==NULL)
		err("cannot open ttfile=%s\n",ttfile);

	if( !getparstring("jpfile",&jpfile)) {
		jpfp = stderr;
	} else  
		jpfp = fopen(jpfile,"w");


	/* get information from the first header */
	if (!fgettr(infp,&tr)) err("can't get first trace");
	nt = tr.ns;
	if (!getparfloat("dt",&dt)) dt = tr.dt/1000000.0; 
	if (dt<0.0000001) err("dt must be positive!\n");
	if (!getparfloat("ft",&ft)) ft = (float)tr.delrt/1000.; 
	
	/* get traveltime table parameters	*/
	if (!getparint("nxt",&nxt)) err("must specify nxt!\n");
	if (!getparfloat("fxt",&fxt)) err("must specify fxt!\n");
	if (!getparfloat("dxt",&dxt)) err("must specify dxt!\n");
	if (!getparint("nzt",&nzt)) err("must specify nzt!\n");
	if (!getparfloat("fzt",&fzt)) err("must specify fzt!\n");
	if (!getparfloat("dzt",&dzt)) err("must specify dzt!\n");
	if (!getparint("ns",&ns)) err("must specify ns!\n");
	if (!getparfloat("fs",&fs)) err("must specify fs!\n");
	if (!getparfloat("ds",&ds)) err("must specify ds!\n");
	ext = fxt+(nxt-1)*dxt;
	ezt = fzt+(nzt-1)*dzt;
	es = fs+(ns-1)*ds;

	if (!getparint("nxso",&nxso)) err("must specify nxso!\n");
	if (!getparfloat("fxso",&fxso)) fxso = 0;
	if (!getparfloat("dxso",&dxso)) err("must specify dxso!\n"); 
        if (!getparint("nxgo",&nxgo)) err("must specify nxgo!\n"); 
        if (!getparfloat("fxgo",&fxgo)) fxgo = 0;
        if (!getparfloat("dxgo",&dxgo)) err("must specify dxgo!\n"); 

        /* errors for illegal survey geometry */        
        if (nxgo>1 && dxgo != dxso) err("in this implementation, dxso must = dxgo!\n");
        if (dxgo < 0.) err("in this implementation, dxgo must be >0!\n");
        if (dxso < 0.) err("in this implementation, dxso must be >0!\n");

        /* calculate the number of independent locations spanning the survey. */
        dxi = dxso;
        fxi = fxso;
        nxi = (int)((fxgo + (nxgo+nxso-2)*dxgo)/dxi) + 1;

        /* calculate the number of independent receiver locations or number of  */
        /* receiver gathers.                                                    */
        fgl = fxso + fxgo;
        ngl = nxgo + (nxso-1); 
 
        /* location of last source and receiver */
        exgo = fxgo + (nxi-1)*dxgo;
        exso = fxso+(nxso-1)*dxso;

        /* get the depth resolution for interploation */
	if (!getparfloat("dzo",&dzo)) dzo = dzt*0.2;
        fzo=fzt;
        nzo=1+(((nzt-1)*dzt)/dzo);
        
        fzo = fzo*sgn;
        dzo = dzo*sgn;
	ezo = fzo+(nzo-1)*dzo;

        /* error if survey dimensions exceed those of the travel time table */
	if(fxt>fxso||fxt>fxgo||ext<exso||ext<exgo||fzt>fzo||ezt<ezo) 
		err(" output range is out of traveltime table!\n");

        if (!getparint("recsurf",&recsurf)) recsurf=0;
        if (!getparint("datsurf",&datsurf)) datsurf=0;

        if(recsurf==0){     
           if (!getparfloat("zrec",&zrec)) err("must specify zrec when recsurf=0!\n");
        }     

        if(datsurf==0){     
           if (!getparfloat("zdat",&zdat)) err("must specify zdat when datsurf=0!\n");
        }    

	if (!getparfloat("v0",&v0)) v0 = 1500;
	if (!getparfloat("aperx",&aperx)) aperx = 0.5*nxt*dxt;
	if (!getparfloat("offmax",&offmax)) offmax = 99999.;

        if (!getparint("dc",&dc))   dc = 0;
        if (dc > 1) err("dc must be 0 or 1!\n");

        if (!getparfloat("freq",&freq))    freq=50;
        if (freq <= 0.0) err("freq must be positive and non-zero!\n");
        if (1/(2*freq)<dt) warn("freq set too high - you may have numerical singularities !\n");

        if (!getparfloat("scale",&scale))    scale=1.0;

	if (!getparint("ntr",&ntr))	ntr = 100000;
	if (!getparint("mtr",&mtr))	mtr = 100;

	fprintf(jpfp,"\n");
	fprintf(jpfp," Datuming parameters\n");
	fprintf(jpfp," ================\n");
	fprintf(jpfp," datain=%s \n",datain);
	fprintf(jpfp," dataout=%s \n",dataout);
	fprintf(jpfp," ttfile=%s \n",ttfile);
        fprintf(jpfp," \n");
        fprintf(jpfp," nzt=%d fzt=%g dzt=%g\n",nzt,fzt,dzt);
	fprintf(jpfp," nxt=%d fxt=%g dxt=%g\n",nxt,fxt,dxt);
 	fprintf(jpfp," ns=%d fs=%g ds=%g\n",ns,fs,ds);
        fprintf(jpfp," \n");
        fprintf(jpfp," ngl=%d fgl=%g nxi=%d fxi=%g dxi=%g\n",ngl,fgl,nxi,fxi,dxi);
        fprintf(jpfp," \n");
	fprintf(jpfp," nzo=%d fzo=%g dzo=%g\n",nzo,fzo,dzo);
	fprintf(jpfp," nxso=%d fxso=%g dxso=%g\n",nxso,fxso,dxso);
        fprintf(jpfp," nxgo=%d fxgo=%g dxgo=%g\n",nxgo,fxgo,dxgo);
        fprintf(jpfp," \n");        
        fprintf(jpfp," nt=%d ft=%g dt=%g \n",nt,ft,dt);
        fprintf(jpfp," freq=%g v0=%g dc=%d\n",freq,v0,dc);
        fprintf(jpfp," aperx=%g offmax=%g\n",aperx,offmax);
        fprintf(jpfp," ntr=%d mtr=%d \n",ntr,mtr);
        fprintf(jpfp," ================\n");
        fflush(jpfp);


        /* read and create recording and datum surfaces */

        if(recsurf==1){
           if(!getparstring("recfile",&recfile))
                err("you chose recsurf=1, so you must specify a recfile!\n");
           recfp = fopen(recfile,"r");
        }

        if(datsurf==1){
           if(!getparstring("datfile",&datfile))
                err("you chose datsurf=1, so you must specify a datfile!\n");
           datfp = fopen(datfile,"r");
        }     

        checkpars();

        szif = ealloc2float(2,nxi);    


        if(recsurf==0){     
          for(i=0;i<nxi;i++) szif[i][0]=zrec;
        }else{

          for(i=0;i<nxi;i++){
           if(fscanf(recfp,"%f",&r1)<0) err("not enough points in recfile!\n");
           szif[i][0]=r1;
          } 
          fclose(recfp);

        }     

        if(datsurf==0){     
          for(i=0;i<nxi;i++) szif[i][1]=zdat;
        }else{

          for(i=0;i<nxi;i++){
           if(fscanf(datfp,"%f",&r2)<0) err("not enough points in datfile!\n");
           szif[i][1]=r2;
          } 
          fclose(datfp);

        }      


        /* allocate array memory */
        ttab = ealloc3float(nzt,nxt,ns);
        tsum = alloc2float(nzt,nxt);
  
        /* array with length equal to number of receiver gathers */ 
        ng = ealloc1int(ngl);
        memset((void *) ng,(int) '\0',ngl*sizeof(int)); 


        /* allocate output array */
        datg = ealloc3float(nt,nxso,ngl);
        memset((void *) datg[0][0],(int) '\0',ngl*nxso*nt*sizeof
           (float));           


 	fprintf(jpfp," input traveltime tables \n");
                       
	/* loop over surface locations and read the corresponding slice  */
        /* of the travel time array.                                     */
	for(is=0; is<ns; ++is){
		nseek = nxt*nzt*is;
		fseek(ttfp,nseek*((long) sizeof(float)),0);
		fread(ttab[is][0],sizeof(float),nxt*nzt,ttfp);
        }



/*  begin main loop.   */

/*  program flow and loop structure:  the outermost (main) loop cycles once   */
/*  through the data traces.  for each trace, loop over the output section,   */
/*  (over an aperture centered on the offset of the data trace), adding to    */
/*  the sum at each point in the output section [ix,it] the value on the data */
/*  trace at the retarded or advanced time (depending on extrapolation        */
/*  direction.  Schematically:                                                */
/*                                                                            */ 
/*      loop over data traces [X]                                             */
/*                                                                            */
/*       loop over output traces [ix] (in aperture centered at X )            */ 
/*                                                                            */
/*         loop over times along output trace [it]                            */
/*                                                                            */
/*            add to sum at output point [ix,it] the value                    */
/*            on data trace at [X,it+-tau]                                    */ 
/*                                                                            */
/*            (where tau is the travel time between the                       */
/*             location [ix,zdatum] on the datum and the trace                */
/*             location )                                                     */
/*                                                                            */


	fprintf(jpfp," start source datuming... \n");
	fprintf(jpfp," \n");
	fflush(jpfp);
  

	jtr = 1;
	ktr = 0;
	fxin = fxso;
	io = 0;
	s = 0.;

        for(is=0; is<ngl; ++is)
                ng[is]=0;               

	/* Store headers in tmpfile while getting a count */
	hdrfp = etmpfile(); 

	do {
     
	    sx = tr.sx;
	    gx = tr.gx;

            /* determine whehter this is the start of a new gather */
	    if(gx!=s){
		fxin = sx;
		s = gx;
	    }

            /* determine which output gather this trace belongs to */
	    io = (int)((gx-fgl)/dxgo);  


	    if(MIN(sx,gx)>=fs && MAX(sx,gx)<=es && 
	        MAX(gx-sx,sx-gx)<=offmax ){

		/* count the number of traces in receiver's gathers  */ 
                ng[io]++;

                efwrite(&tr, 1, HDRBYTES, hdrfp);

                /* interpolate travel time table to the source location */
	    	as = (sx-fs)/ds;
	    	is = (int)as;
		if(is==ns-1) is=ns-2;
		res = as-is;
		if(res<=0.01) res = 0.0;
		if(res>=0.99) res = 1.0;
		sum2(nxt,nzt,1-res,res,ttab[is],ttab[is+1],tsum);


                /* call to dat2d will take the energy on this data trace */
                /* and remap it to the output along the appropriate     */
                /* operators, after taking the required derivatives and */
                /* multiplying by the appropriate scale factors.        */


                dat2d(tr.data,nt,ft,dt,sx,gx,datg[io],aperx,nxso,fxin,
                      dxso,nzo,fzo,dzo,nxi,fxi,dxi,
                      tsum,sgn,nzt,fzt,dzt,nxt,fxt,dxt,
                      szif,v0,dc,ngl,freq);


	        ktr++;
	        if((jtr-1)%mtr ==0 ){
			fprintf(jpfp," Datumed source in trace %d\n",jtr);
			fflush(jpfp);
	    	}

	    }
	    jtr++;
	} while (fgettr(infp,&tr) && jtr<=ntr);

/*  end of main loop  */


/*  output data  */

	fprintf(jpfp," Datumed sources in %d total traces\n",ktr);

	/* Seg-Y header	*/

	rewind(hdrfp);


	/* scale factors representing the constant in the integral */

        scal = dxso*(1/sqrt(2*PI*v0));
        scal = scal*scale;


        for(io=0; io<ngl; io++) {
        	for(ixo=0; ixo<ng[io]; ixo++)  {

			efread(&tro, 1, HDRBYTES, hdrfp);
			memcpy((void *) tro.data,
			   (const void *) datg[io][ixo], nt*sizeof(float));

			for(it=0; it<nt; it++){
			   tro.data[it] *=scal;
                        }

			/* write out */
			fputtr(outfp,&tro); 

		}
	}

	fprintf(jpfp," \n");
	fprintf(jpfp," output done\n");
	fflush(jpfp);

	efclose(jpfp);
	efclose(outfp);
	efclose(hdrfp);

	    
        free2float(tsum);
        free3float(ttab);
        free3float(datg);
	return EXIT_SUCCESS;
}

/* sum of two tables	*/
  void sum2(int nx,int nz,float a1,float a2,float **t1,float **t2,float **t)
{
	int ix,iz;

	for(ix=0; ix<nx; ++ix) 
		for(iz=0; iz<nz; ++iz)
			t[ix][iz] = a1*t1[ix][iz]+a2*t2[ix][iz];
}
 

  void dat2d(float *trace,int nt,float ft,float dt,float sx,float gx,
             float **dat,float aperx,int nxso,float fxin,float dxso,
             float nzo,float fzo,float dzo,
             int nxi,float fxi,float dxi,
	     float **tsum,int sgn,
             int nzt,float fzt,float dzt,int nxt,float fxt,float dxt,
             float **szif,float v0,int dc,int ngl,float freq)

/*****************************************************************************
Datum source on one trace 
******************************************************************************
Input:
*trace		one seismic trace 
nt		number of time samples in seismic trace
ft		first time sample of seismic trace
dt		time sampleing interval in seismic trace
sx,gx		lateral coordinates of source and geophone (ignored)
aperx		lateral aperature in datuming
nxso,fxin,dxso  traces in output 
nzo,fzo,dzo     interpolation variables to define surfaces	
fxi             x-coordinate of the first surface location
dxi             horizontal spacing on surface
nxi             number of surface locations
tsum            sum of residual traveltimes from shot and receiver
nxt,fxt,dxt,nzt,fzt,dzt        dimension parameters of traveltime table
szif            array[] of the z component of the reference surfaces
dc              flag for previously downward continued receivers (dc=1)
ngl             number of geophone locations = number of receiver gathers

Output:
dat		Redatumed section
*****************************************************************************/
{
	int nxf,nxe,it,ix,iz,jz,jt,jx,mr;
	float x,z,ampd,res0,am,am0,ax,ax0;
	float odt=1.0/dt,az,sz,sz0,at,tio,res,temp,fxf,h;
        float p,q,s,ctau,xst=0.0,zst=0.0,ris,ros,rg,t0=0.0,sing;
        float xp,zgx,hp,phi=0.0,sxp,zs,zi,zip,atol;
	float *tzt;
        float zp=0.0,dzde=0.0,g,xs;


        /* set a tolerance variable to deal with numerical problems */
        atol = 1e-10; 

	/* allocate memory for datuming */
	tzt = alloc1float(nzt);

        h = (gx>=sx)?(gx-sx)/2.0:(sx-gx)/2.0;
	fxf = fxi + (nxi-2)*dxi;

        /* call to filt will perform the frequency domain filtering  */
        /* on the input data trace, returning in variable trace.     */  

        filt(trace,nt,dt,sgn);                


        /* calculate the limits of the aperture in the output */
	nxf = (sx-aperx-fxin)/dxso;
	if(nxf<0) nxf = 0;
        nxe = (sx+aperx-fxin)/dxso;
	if(nxe>=nxso) nxe = nxso-1;


        /* error if survey length exceeds topography definition */ 
        am = (fxin-fxi)/dxi;
        mr = (int)am;
        am = am - mr;
        if(am<=0.01) am = 0.;
        if(am>=0.99) am = 1.0;
        am0 = 1.0-am;
        if(mr<0) mr = 0;
        if(mr+nxe>nxi-1)
                err("Topography definition is out of range!\n");
	

        /* depth of the input source(zi) on the recording surface */

        am = (sx-fxi)/dxi;
        mr = (int)am;
        am = am - mr;
        if(am<=0.01) am = 0.;
        if(am>=0.99) am = 1.0;
        am0 = 1.0-am;
        if(mr<0) mr = 0;

        if(mr>=nxi-1){
           mr=mr-1;
           am0 = 0.0;
           am = 1.0;
        }

        zi = ABS(am0*szif[mr][0]+am*szif[mr+1][0]);
        if(zi>=(nzo*dzo)+fzo)
              err("Recording surface is out of range!\n");


        /* approximate the surface derivative near this input location */

        if(mr==0){ 
           zp = szif[mr+1][0]-szif[mr][0];
           dzde = zp/dxi;         
        }else{
           zp = (szif[mr][0]-szif[mr-1][0]) + (szif[mr+1][0]-szif[mr][0]);
           dzde = zp/(2*dxi);
        }


        /* depth of the geophone location(zgx), depending on whether */
        /* it has previously been downward continued.  if not, put it   */
        /* on the recording surface, if so, put it on the datuming surface. */

        am = (gx-fxi)/dxi;
        mr = (int)am;
        am = am - mr;
        if(am<=0.01) am = 0.;
        if(am>=0.99) am = 1.0;
        am0 = 1.0-am;
        if(mr<0) mr = 0;

        if(mr>=nxi-1){
           mr=mr-1;
           am0 = 0.0;
           am = 1.0;
        }

        if(!dc){

            zgx = ABS(am0*szif[mr][0]+am*szif[mr+1][0]);
            if(zgx>=(nzo*dzo)+fzo)
                  err("Recording surface is out of range!\n");

        }else{

            zgx = ABS(am0*szif[mr][1]+am*szif[mr+1][1]);
            if(zgx>=(nzo*dzo)+fzo)
                  err("Datuming surface is out of range!\n");
        }


        /* loop over the lateral aperture in the output section */
	for(ix=nxf; ix<=nxe; ++ix){
		x = fxin+ix*dxso;

		ax = (x-fxt)/dxt;
		jx = (int)ax;
		ax = ax-jx;
		if(ax<=0.01) ax = 0.;
		if(ax>=0.99) ax = 1.0;
		ax0 = 1.0-ax;


                /* determine the time tio from the source location of this */
                /* input trace to the source location of this output trace */
                /* by building a vector of times to all depths, then interp- */
                /* -olating to the exact datuming depth.                     */

		for(iz=t0; iz<nzt; ++iz){
		    tzt[iz] = ax0*tsum[jx][iz]+ax*tsum[jx+1][iz];
		}


                /* datum depth at this output source location(z) */
  
                am = (x-fxi)/dxi;
                mr = (int)am;
                am = am - mr;
                if(am<=0.01) am = 0.;
                if(am>=0.99) am = 1.0;
                am0 = 1.0-am;
                if(mr<0) mr = 0;

                if(mr>=nxi-1){
                     mr=mr-1;
                     am0 = 0.0;
                     am = 1.0;
                }
 
                z = ABS(am0*szif[mr][1]+am*szif[mr+1][1]);
                if(z<=fzo || z>=fzo+(nzo*dzo))                                                    
                      err("Datuming surface is out of traveltime range!\n");    



                /* calculate travel time between datum and surface */
                az = (z-fzo)/dzt;
                jz = (int)az;
                if(jz>=nzt-1) jz = nzt-2;
                sz = az-jz;
                sz0 = 1.0-sz;
                tio = sz0*tzt[jz]+sz*tzt[jz+1];


/* add this for exact straight ray calc */

/*
                tio =  sqrt((x-sx)*(x-sx)+(z-zi)*(z-zi))/v0;
*/

/*                                      */


                at = (sgn*tio-ft)*odt;

                /* loop over times along the output trace, and add to the */
                /* sum at the output point [ix,it] any amplitude on the   */
                /* data trace at time it+-tio.                            */

		jt = (int)at;
                res = ABS(at-jt);
                res0 = 1.0-res;


                /* only sources on the lesser side of this receiver location are in */
                /* this gather by the required geometry.  if the output location is */
                /* not, do not use it.                                              */
 
                if(x<=gx){

                  /* define the time along the singular path */  
                  sing = sqrt((gx-x)*(gx-x)+(z-zgx)*(z-zgx))/v0;

                  /* define the rotation */
                  if(x==gx){
                    if(dc==0) phi = PI/2.0;
                    if(dc==1) phi = 0.0;
                  }else{
                    phi = atan((z-zgx)/(x-gx));      
                  }

                  /* coordinate shift-rotation */
                  xp = (x-gx)*cos(phi)+(z-zgx)*sin(phi);
                  hp = 0.5*xp;
                  sxp = (sx-gx)*cos(phi)+(zi-zgx)*sin(phi);
                  zip = -(sx-gx)*sin(phi)+(zi-zgx)*cos(phi);

                  if(dc==0 && hp==0.0) err("Datum and recording surface overlap!\n");


                  /* loop over output times */
                  for (it=0; it<nt; it++){
		    if(it+jt >= 0 && jt < nt-it-sgn){
  
                     t0 = ft+(it*dt); 

                     if(t0<sing+atol){ 

                        ampd=0.0;

                     }else{                
  
                        /* calculate the coordinates of the stationary point */
                        /* in the shifted-rotated frame */


                        if(dc==0){

                           if(h==0.0){

                              xst=hp+v0*t0/2.0; 
                              zst=0.0;

                           }else{

                              ctau = (v0*t0)*(v0*t0); 
                              q = (ctau-4.0*hp*hp)/(4.0*ctau); 

                              p = 2.0*q*(sxp-xp)*(xp)*zip;      
                              p = p/((zip*zip) + 4.0*q*(sxp-xp)*(sxp-xp));    

                              s = q*zip*zip*(ctau - (xp*xp));          
                              s = s/((zip*zip) + 4.0*q*(sxp-xp)*(sxp-xp));                           

                              zst = sqrt(s+(p*p))-p;        
                              xst = ((sxp-xp)/zip)*zst + xp;

                              if(zst<0.0) err("Negative zst!\n"); 

                           }
                        }

                        if(dc==1 && hp!=0.0){

                           if(h==0.0){

                              xst=hp+v0*t0/2.0; 
                              zst=0.0;

                           }else{

                              ctau = (v0*t0)*(v0*t0); 
                              q = (ctau-4.0*hp*hp)/(4.0*ctau); 

                              p = 2.0*q*(sxp-xp)*(xp)*zip;      
                              p = p/((zip*zip) + 4.0*q*(sxp-xp)*(sxp-xp));    

                              s = q*zip*zip*(ctau - (xp*xp));          
                              s = s/((zip*zip) + 4.0*q*(sxp-xp)*(sxp-xp));                  

                              zst = sqrt(s+(p*p))-p;        
                              xst = ((sxp-xp)/zip)*zst + xp;

                              if(zst<0.0) err("Negative zst!\n"); 

                           }
                        }
                            
                        if(dc==1 && hp==0.0){

                              xst = -(0.5)*(t0*v0*sxp)/sqrt((sxp*sxp)+(zip*zip));  
                              zst = ABS((0.5)*(t0*v0*zip)/sqrt((sxp*sxp)+(zip*zip)));        

                        }


                        /* rotate back to unprimed coords */

                        xs = xst*cos(phi)-zst*sin(phi)+sx;
                        zs = xst*sin(phi)+zst*cos(phi)+zgx;           
                  /*    if(zs<0) err("Negative zs...\n");    */

                        /* calculate path lenghts */

                        ris = sqrt((xst-sxp)*(xst-sxp) + (zst-zip)*(zst-zip));
                        ros = sqrt((xst-xp)*(xst-xp) + (zst)*(zst));     
                        rg = sqrt((xst)*(xst) + (zst*zst));    

                   /*   if(ris<=0. || ros<=0. || rg<0.) warn("Negative path length!\n");  */
                   /*   if(ris<=ros) warn("St. pnt. not below output!\n");   */

                        /* calculate G(xs,zs) */ 

                        g = (zs-zi)-dzde*(xs-sx);

                        /* calculate the amplitude weight */

                        ampd = g*sqrt(rg+ris)/ris;                 
                        ampd = ampd/(sqrt(rg+ros)*sqrt(ABS(ris-ros)));                       

                        /* check to see if this is a valid stationary point */
                        if((ris-ros)<(v0/(4*freq))) ampd = 0.0;                           
                        if(zs<0.) ampd = 0.0;
                        if(ris<=ros) ampd=0.0;
                        if(ris<0. || ros<0. || rg<0.) ampd=0.0;

                     }             

  
                     /* find the correct time data, scale by amplitude, and */
                     /* assign to output array. */ 

                     temp=(res0*trace[it+jt]+res*trace[it+jt+sgn])*ampd;
                     dat[ix][it] += temp;       


                    }
              	}

              }
	}

	free1float(tzt);
}

void filt(float *trace,int nt,float dt,int sgn)

/*******************************************************************
  input: 
    trace(nt)   single seismic trace
  output:
    trace(nt)   filtered and phase-shifted seismic trace 
********************************************************************/
{
        static int nfft, n;
        static float dw, fw;
        int  it,iw,sgnw;
        float temp, const2, amp, omega;
        complex *ct;


        /* later calculations will require multiplication by sqrt(-i) */
        /* which is equivalent to a phase shift of pi/4.  it will be  */
        /* applied in Euler form (cos-isin), so define a constant:    */
        /*                                  const2=cos(pi/4)=sin(pi/4)*/
        const2 = 0.5*sqrt(2.0);

        if(nfft==0) {   
                /* Set up FFT parameters */

                nfft = npfao(2*nt,4*nt);              

                if (nfft >= SU_NFLTS || nfft >= 720720)
                        err("Padded nt=%d -- too big", nfft);

                n = nfft/2;
                dw = 2.0*PI/(nfft*dt);
                fw = -PI/dt;

        }        

        /* allocate complex array for fft */

        ct   = ealloc1complex(nfft);

        /* zero the real and complex parts of ct out to nfft */
        for(it=0; it<nfft; ++it){
             ct[it].r = 0.0;
             ct[it].i = 0.0;
        }
 
        /* put the trace values into the real part of ct- already padded */
        for(it=0; it<nt; ++it){
             ct[it].r = trace[it];
        }


        /* do the transform for t to omega */
        pfacc(1,nfft,ct);



        /* do omega-domain filtering for point source - 2.5D */

           /* do positive frequencies */
           for(iw=0; iw<=n; ++iw){                           
               
                omega = iw*dw; 
                amp = sqrt(ABS(omega));
                sgnw = 1;    

                temp = (ct[iw].r-sgnw*ct[iw].i)*amp*const2;        
                ct[iw].i = (sgnw*ct[iw].r + ct[iw].i)*amp*const2;  
                ct[iw].r = temp;                                      
           }                                 

           /* do negative frequencies */
           for(iw=n+1; iw<=nfft-1; ++iw){

                omega = fw+(iw-n)*dw;
                amp = sqrt(ABS(omega));
                sgnw = -1;     

                temp = (ct[iw].r-sgnw*ct[iw].i)*amp*const2;        
                ct[iw].i = (sgnw*ct[iw].r + ct[iw].i)*amp*const2;  
                ct[iw].r = temp;                                      
           }                               


        /* do the inverse transform for omega to t */

        pfacc(-1, nfft, ct);            
                
        for (it=0; it<nt; ++it) {
             trace[it] = ct[it].r/nfft;
        }

        free1complex(ct);
}

